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Unresolved near-infrared background anisotropies are expected to have contributions from 
the earliest galaxies during reionizationP^and faint, dwarf galaxies at intermediate redshifts^l 
Previous measurements™^ were unable to conclusively pinpoint the dominant origin because 
they did not sample spatial scales that were sufficiently large to distinguish between these 
two possibilities. Here we report a measurement of the anisotropy power spectrum from 
sub-arcminute to one degree angular scales and find the clustering amplitude to be larger 
than the model predictions involving the two existing explanations. As the shot-noise level of 
the power spectrum is consistent with that expected from faint galaxies, a new source pop- 
ulation on the sky is not necessary to explain the observations. A physical mechanism that 
increases the clustering amplitude, however, is needed. Motivated by recent results related to 
the extended stellar light profile in dark matter halos^SSS we consider the possibility that the 
fluctuations originate from diffuse intrahalo stars of all galaxies. We find that the measured 
power spectrum can be explained by an intrahalo light fraction of 0.07 to 0.2% relative to 
the total luminosity in dark matter halos of 10 9 to 10 12 solar masses at redshifts of ~ 1 to 4. 

In order to distinguish between the two interpretations of the near-IR anisotropy power spec- 
trum, we have analyzed imaging data from the Spitzer Deep, Wide-Field Survey (SDWFS} 16 . This 
survey covers 10.5 square degrees on the sky with the IRAC instrument in its four bands between 
3.6 and 8 /im. We focus on the data at 3.6 and 4.5 microns as the confusion from zodiacal light 
limits extragalactic background studies at 5 and 8 microns®. The data were taken in four separate 
epochs between 2004 and 2008 and were conducted in ways to minimize the systematics associated 
with anisotropy measurements. In particular the four different epochs were observed at different 
roll angles of the instrument so that the measurements are robust against detector artifacts, persis- 
tence resulting from saturated bright stars, and variations in the bias level. The SDWFS mapping 
strategy was also optimized to facilitate self-calibration^ of the data by maximizing inter-pixel 
correlations. 
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To limit the influence of bright stars and galaxies, including extended sources, in our anisotropy 
measurements we mask all sources that are detected either in the combined SDWFS data or in the 
ancillary multi-band optical and near-IR data 18 . The effects of mosaicing of individual detector 
frames, pixelization of the maps, and the detected-source mask are captured by the map-making 
transfer function (see Supplementary Information). We compute the transfer function and its un- 
certainty with a large set of sky simulations. The point-spread function (PSF) and its uncertainty 
were determined by measuring and modeling the PSF of stars at different sub-regions of the image 
and computing the variance of the differences between the modeled PSFs. 

The power spectrum measurements at 3.6 and 4.5 fxm show a clear excess above the shot- 
noise level (Fig. 1). The shot-noise dominates the anisotropy power spectrum at sub-arcminute 
angular scales corresponding to £ > 10 5 . Such a shot-noise is expected from the small-scale Pois- 
son behavior of the spatial distribution of sources on the sky. The clustering amplitude we measure 
at £ ~ 10 4 is fully consistent with existing measurements of the anisotropy power spectrum with 
IRAC in deeper, but smaller area, fields^^l. At tens of arcseconds angular scales, corresponding 
to £ > 5 x 10 4 , our shot-noise level is higher than that of a recent measurement by about a factor 
of 2 because deeper data allow more faint sources to be individually detected and masked 12 . Nev- 
ertheless, we independently confirm the near-IR background anisotropics at angular scales larger 
than a few arcminutes at the previously reported amplitude^. 

The near-IR anisotropics have been previously interpreted as either due to spatial clustering 
of primordial galaxies responsible for cosmic reionization 8 or due to faint, dwarf galaxies at low 
redshifts that fall below the individual source detection threshold of Spitzer images^^l Fig. 1 
shows that the measured fluctuations in SDWFS are well above both these interpretations. The 
power spectrum predictions for z > 6 galaxies rely on a combination of analytical calculation^ 
and numerical simulations 20 of reionization. If we force the z > 6 galaxy model to fit the power 
spectrum data then the integrated intensity of z > 6 galaxies is about 2 nW m~ 2 sr^ 1 at 3.6 /xmPl 
In order to reach such a high intensity these galaxies must be very efficient in converting baryons 
to starP^. In fact, the required star-formation rate conflicts with the measured metal abundance at 
z > 4, the measured X-ray background when compared to X-rays from stellar end products such 
as black holes, and the measured luminosity functions of bright Lyman-dropout galaxies^. Unless 
a significant revision of our current understanding of z > 6 galaxy statistics is made it is unlikely 
that the measured anisotropy power spectrum is dominated by the primordial galaxies. 

The prediction for low redshift, faint galaxy intensity fluctuations involves a large compila- 
tion of multi-wavelength luminosity functions and galaxy number counts^. The measured luminos- 
ity function slope at the faint-end is used to extrapolate to the fainter galaxies that are undetected in 
the Spitzer images. An increase in the faint-end slope above the measured values does not increase 
the clustering amplitude on the angular scales of interest without modifying the shot-noise level. 
While the clustering amplitude is smaller than the measurements, the prediction related to faint, 
dwarf galaxies 7 shows that they generate a shot-noise level consistent with the measured small- 
scale anisotropy power spectrum (Fig. 1). At few tens arcminute angular scales the measurements 
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are such that the clustering amplitude is about a factor of 6 to 10 above the prediction. While this 
difference suggests that a new model to explain the anisotropy power spectrum is clearly needed, 
the consistency with the shot-noise level is such that we do not need to invoke a new population of 
point sources on the sky to explain the observations. 

While keeping the shot-noise level the same, the measurements can be explained by any 
physical effect that boosts the two-halo term of clustering. One possibility is to increase the halo 
mass scale of the faint, dwarf galaxies so that clustering bias factor is increased. The required 
modification needed to explain the fluctuations data, however, is ruled out by the measured number 
counts and the redshift distribution 7 . Since intensity anisotropics are measured, another option is to 
introduce a luminosity component to the dark matter halos that remain unmasked when the hosted 
bright galactic disks are masked as part of the analysis. Such a possibility exists in the literature in 
the form of diffuse halo stars in the extended stellar profile of galaxies out to distances of 100 kpcPl 
In our anisotropy measurements, we mask the faintest detected galaxies to 3-4" which removes the 
light from the bulges and disks of those galaxies. To remove the diffuse light component we would 
have to mask to a radius greater than 10" around each galaxy. The surface density of galaxies down 
to ttiab < 22 at 3.6 /mi is such that we expect 2 to 3 galaxies within a circle of radius 10". Thus 
masks which successfully remove the diffuse component leave no pixels on the map from which 
to measure the anisotropy power spectrum. 

Existing studies discuss this extended emission in terms of the diffuse intrahalo light (IHLj^ 
and explain the origin through tidally stripped stars during galaxy mergers and collisions. The 
stripped fraction is expected to be a function of the halo mass with more massive halos containing 
a larger fraction of the diffuse halo emissiorP^SEl On galaxy cluster scales the diffuse intracluster 
light ! 24 * 25 ! is a significant fraction of the total luminosity of the cluster. We describe the intensity 
anisotropy power spectrum from the IHL by modifying the standard galaxy clustering predictions^ 
to include a profile for the diffuse stars in halos (see Section 8 of the Supplementary Information). 
If the clustering excess in near-IR anisotropy power spectrum is due to IHL, then we find that 
measured anisotropics can be described with halos in the mass range of 10 9 to 10 12 M . Averaged 
over this mass range we find an fmh of 0.07 to 0.2% at 68% confidence level (Fig. 2). The implied 
fraction is consistent with the theoretical expectation that the IHL level is small for low mass halos, 
but differences also exist with current theory prediction s ^ 13 * 14 !, especially in terms of the power-law 
slope of the halo mass dependence. 

If this new interpretation involving IHL is the correct description of measured IR background 
anisotropics, we find that the IHL in all dark matter halos that we are probing contribute 0.75±0.25 
nW m~ 2 sr _1 to the total intensity at 3.6 yum. This intensity can be compared to the rms fluctuations 
of about 0.1 nW m~ 2 sr _1 at a few arcminutes angular scale (see Fig. 3). The IHL fluctuation 
signal varies spatially at a level of 10 to 15% of its integrated intensity and below 1% of the total 
background intensity of 13.3 ± 2.8 nW m~ 2 sr _1 at 3.6 /imPl As the spectral energy distribution 
of IHL is mostly unknown, we make use of a variety of SEDs from B to K-type stellar spectral 
templates and find an order of magnitude variation at wavelengths of ~ 1 /im (Fig. 3). In all 
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these cases we predict the existence of optical background light fluctuations. They will have a 
similar power spectrum shape and will be fully correlated with fluctuations at 3.6 /im. Furthermore 
the near-IR anisotropics we have measured should be correlated with the sub-mm anisotropics^, 
especially if there are diffuse and extended dust associated with galaxies. These form future tests 
that can be used to improve our understanding of the content and nature of IHL in distant dark 
matter halos. 
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Figure 1: The angular power spectrum of the unresolved near-IR background. The total 
power spectrum P(k) of SDWFS at 3.6 /im (top) and 4.5 /im (bottom) as a function of the mul- 
tipole moment. SDWFS imaging data were taken on the same field at four separate epochs in 
January 2004, August 2007, February 2008, and March 2008. Each epoch of data, taken over 7 
to 10 days, includes 4300 to 4900 IRAC frames that were combined to make mosaics using the 
self-calibration algorithm 17 . The total integration time is 6 minutes per pixel. These individual 
frames were first visually inspected and cleaned of artifacts such as asteroidal trails and hot pixels. 
Through cross-correlations between sum and difference maps between epochs, we make indepen- 
dent measurements of the sky signal and noise. The final power spectrum is the average of the 
multi-epoch cross-correlation data under the assumption that the instrumental noise is not corre- 
lated between epochs. In both panels the error bars are la uncertainties in the power spectrum. 
They are determined by propagating the errors from the beam measurement into the power spec- 
trum, while the simulations, based out of noise measurements, were used to obtain instrumental 
and sky variance. The quadratic sum of these errors and the map-making transfer function uncer- 
tainty constitutes the final error estimate. The two shaded regions show the expected contribution 
from z > 6 galaxie^ and low-redshift galaxieP based on two model predictions in the literature. 
The lines shows a diffuse intrahalo light model where we show the signal in terms of the total 
(solid), one (dashed-dotted) and two (dotted) halo terms. The dashed line is the best-fit shot-noise 
signal that dominates the anisotropics at small angular scales. 
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Figure 2: The intrahalo light fraction from diffuse stars as a function of the halo mass. The 

dark and light blue shaded regions show the 95% and 68% range of /ihl relative to the total 
luminosity of the dark matter halos as a function of the halo mass from an analytical prediction^ 
valid for f mL > 4 x 10~ 4 and M > 5 x 10 10 M and at z = 0. We show the case where subhalos 
on orbits passing within a critical radius of the host halo center contribute their light to the central 
galaxy rather than to the diffuse component. We also show a prediction where /ihl is constant, 
due to dwarf galaxies that are completely destroyed, with a value ~ 0.005 when M < 5 x 10 11 
M Q (solid line fixed at / ffl L = 5 x 10~ 3 ). The downward arrow indicates the possibility that the 
constant /ihl value for low mass halos may be smaller at higher redshifts. The red and orange 
hatched regions at the bottom of the plot are the preferred 68% and 95% confidence level range, 
respectively, on /ihl from our analysis of the SDWFS near-IR anisotropy power spectrum. The 
mass range is determined by the minimum and maximum halos masses consistent with the halo 
model fit that includes the IHL component. Both the mass and /ihl ranges are valid over the broad 
redshift interval from z = 1 to 4 over which the anisotropy signal is generated. We do not find 
a significant halo mass dependence on the IHL fraction with the mass -dependent power-law to be 
0.09 ±0.01 between 10 9 to 10 12 M Q , consistent with the possibility that /ihl is mass independent!^ 
when M < 5 x 10 11 M . Our model requires the total luminosity-halo mass relation to evolve with 
redshift as (1 + z) L2±al . This luminosity evolution with redshift can also be absorbed into the 
evolution of fmh(M) evolution with redshift. For reference, we also show measurements and la 
errors of the intracluster light 23 , the galaxy group and cluster analog for IHL when M > 5 x 10 13 
M Q . At halo masses around 10 12 M we show the 95% confidence level upper limit on /i H l 
estimated for Milky Way^and Andromeda (M31)P. 
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Figure 3: The spectral energy distribution of IR background anisotropics. The frequency 
spectrum of near-IR and optical background anisotropics as a function of the wavelength. We 
show the rms fluctuation amplitude at t = 3000, corresponding to fluctuations at 6 arcminute 
angular scales. We show our measurement and the existing measurements in the literature at the 
same angular scales^S. The plotted error bars are the la uncertainties directly propagated to a 
rms fluctuation amplitude from the power spectrum errors. We do not show the measurements at 
1.6 and 1.1 /im^as they do not probe fluctuations on scales greater than 2 arcminutes on the sky 
due to the narrow field of view of the observations. The shaded regions are (a) predictions for 
the IHL taking a variety of spectral energy distributions for the stripped stars (green), (b) z > 6 
galaxies (yellow), and (c) low-redshift galaxies below the detection level of masking (blue). 
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Supplementary Information 



In these Supplemental Notes to the main paper, we outline key details related to how we 
estimated the angular power spectrum of Spitzer/IRAC data and its interpretation. The data used 
for this study are publicly available from the Spitzer Heritage ArchiveQ under the program number 
GO 40839 (PI. D. Stern). 

1 The Spitzer Deep Wide Field Survey 

We use the IRAC data from the Spitzer Deep Wide Field Survey (SDWFS^™. Of the NDWFS 
Bootes fields, the observations were obtained in 4 epochs, with depth per pixel of 90 seconds in 
each epoch. Each epoch took observations over 7 to 10 days to complete the full mosaic. We 
started with the basic calibrated data (BCDs). At each epoch, imaging data was obtained in all four 
IRAC wavelengths or channels (3.6, 4.5, 5.8, and 8 /im); however, for this study we focus on the 
3.6 and 4.5 /im data. Each epoch consists of 4300 to 4900 BCDs per channel that were mosaiced 
to form the final images used in the fluctuation analysis. 

Instead of using the public SDWFS mosaics (Available at http://irsa.ipac.caltech.edu), made 
with the standard Spitzer data analysis pipeline MOPEXpl, we produced our own mosaics in order 
to better control systematic errors. To mosaic the BCDs for each channel and epoch we used a 
self-calibration algorithrri^Uto properly match the sky background level from one adjacent frame 
to the other in the overlapping region using an optimized least squares fitting technique. The 
SDWFS mapping strategy incorporates several elements to facilitate self-calibration of the data by 
maximizing inter-pixel correlations. We dithered the observations on small scales and offset by 
one-third of an IRAC field-of-view between successive passes through each group. This provides 
inter-pixel correlation information on both small and large scales, so the self-calibrated mosaic 
has background levels that are stable across the wide area of the SDWFS mosaic. Finally, for the 
larger, rectangular groups, we cadence the observations such that revisits cover the same area but 
with a different step size. With a < 10% penalty in mapping efficiency, cadencing significantly 
enhanced the inter-pixel correlations across all scales. 

These mapping strategies were designed to significantly enhance the self-calibration of the 
data. Finally, by reobserving the field multiple times at different roll angles, our observing strategy 
was designed to be robust against bad rows/columns, large scale cosmetic defects on the array, 
after-images resulting from saturation due to bright stars, variations in the bias level, and the color 
dependence of the IRAC flat-field across the array^. In particular, the challenging diffuse back- 
ground measurements we report here are vastly aided by the redundant coverage: independent data 
sets of the same region are the best way to assess and control systematic errors. With this mapping 
strategy we were able to construct independent sky realizations for carrying out jackknife testing. 

1 http://sha.ipac.caltech.edu/applications/Spitzer/SHA/ 
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Furthermore, the power spectrum is estimated by cross -correlating the maps from different epochs, 
eliminating bias from uncorrelated signals such as instrumental noise and mosaicing artifacts. 

For each IRAC channel and epoch we passed the cleaned BCD into our self-calibration code, 
a slightly modified version of an existing codef^, as inputs the cleaned BCDs (cBCDs) with the 
final, output array size and astrometry defined to correspond to the mosaic of first 3.6 /mi epoch. 
The cBCDs were first cleaned of asteroid trails, hot pixels, and other image artifacts. Since the 
astrometry is the same for each mosaic, they can be properly coadded and jack-knifed as described 
below. Each of the cBCDs have an angular pixel scale of 1.2 arcseconds, which was preserved 
in the final mosaics. The portion of the maps used for analysis are ~ 3.5 x 3 degrees on a side 
for a total area of ~ 10.5 square degrees. The final mosaics generated from the self-calibration 
algorithm are shown in Fig. SQ3 



SDWFS Bootes 3.6 nm (nW/m 2 Sr) SDWFS Bootes 4.5 (nW/m 2 Sr) 




0.5 1.0 1.5 2.0 2.5 3.0 0.5 1.0 1.5 2.0 2.5 3.0 

(Degrees) ; (Degrees) 



Figure S 1: The SDWFS Maps. The final 3.6 /mi (left) and 4.5 /im (right) self-calibrated mosaics 
used for the Bootes SDWFS analysis. 

2 Generation of the Detected Source Mask 

As the analysis concentrates on the background intensity fluctuations, with the aim of identifying 
the nature of faint sources below the individual detection level, we must remove the contamination 
from individual detected sources in our power spectrum measurement. We created source masks 
based both on the objects detected in the IRAC data at 3.6 and 4.5 /im and the NOAO Deep Wide 
Field Survey (NDWFS) catalogs for the B^,, R and I bands. The masked sources include stellar 
point sources, galaxies that are extended, and galaxies that are unresolved but detected as point 
sources in SDWFS. The mask also accounts for the Spitzer-IRAC point spread function (PSF). For 
this study we make use of the publicly available "extended" IRAC PSF to properly account for flux 
wings. 
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Next, we summarize our recipe to generate the source list and discuss more details below: 

1. We create a catalog consisting of all objects detected by Sextractor at a 3a detection 
threshold on the coadded and individual epoch maps at both 3.6 and 4.5 yum. 

2. We add any additional objects detected in the NDWFS optical catalogs, but missing from the 
IRAC catalogs. 

3. We create a map initially of zeros where we place each source detected from our combined 
catalogs with the proper flux values. The sources that are flagged as extended in the SExtrac- 
tor catalogs are placed as extended sources with a size that is comparable to their estimated 
size. 

4. We convolve this final map with the IRAC PSF to ensure that the flux wings for each source 
are properly masked. 

5. We make a cut at a certain flux level so that all pixels with intensity down to that threshold 
are masked. We then histogram the remaining pixels from the data map and cut all pixels that 
are ± 5 a away in the histogram. The latter step is similar to prior approaches to measure 
the Spitzer-IRAC background anisotropy power spectrurrP^. 

6. To produce the final mask, we set the pixels with intensity values above the final flux and out- 
liers from histogram to a value of zero, and the remaining pixels, the ones used for anisotropy 
measurements, a value of one. 

Before the mosaics were created, source extractor was run iteratively on each epoch and 
waveband individually as well as on the coadded maps for both wavebands in order to find detected 
sources. The parameters used for Sextractor are the same as the ones used for the original 
SDWFS catalogs 16 . The combined catalog obtained from this iterative source extractor analysis, 
as well as the objects detected in the NDWFS catalog are merged into a final catalog. 

The final flux cut for the mask was found iteratively by lowering the flux limit until fur- 
ther expansions of the masked regions no longer affect the final power spectrum. Since the PSF 
is slightly different for the 3.6/im and 4.5/xm wavebands the two masks are independently con- 
structed. Fig. S[2] shows the final masked epoch 1 maps. The final mask is such that 56% of the 
original pixels are removed from the subsequent analysis. 

3 Power Spectrum Estimation 

With the final mosaics and masks in hand, initial, raw auto and cross-correlations may be computed 
to measure the level of clustering in the maps at each scale. To calculate the cross -correlation 
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Figure S 2: The Masked SDWFS maps. The final 3.6 /im (left) and 4.5 /im (right) masked maps 
used for the Bootes SDWFS analysis. The mask used for the analysis removes 56% of the pixels 
in the mosaic. 



between masked maps M\ and M 2 in real space, we first take the 2D Fourier transform of each 
map which we call Mi and M 2 , as 

Af-l N-l 

M[l x J y ] = A J2 E M{m,n]e- 2ni{l * m/M+l y n/N) (1) 

m=0 n=0 

where M and N are the number of discrete points in the two dimensions of the map and A is the 
sampling interval in radians. The power spectra C\ formed from the cross-correlation of M\ and 
M 2 for a specific U bin between between /-modes l\ and l 2 is equal to the weighted mean of the 
squared Fourier modes MiM 2 between li and l 2 , 

^ll+ll>h w [^ly\ 

where w[l x , l y ] is a window function in Fourier space that is non-zero for each mode of the analysis 
and zero for modes that are discarded. To compute the raw auto-correlation we have the special 
case where Mj = M 2 . 



Being able to form a power spectrum from a cross-correlation rather than an auto-correlation 
is highly advantageous, as the noise bias and other contaminants that can dominate an auto- 
correlation calculation are minimized in a cross -correlation. This is because the pixel of each 
map Mi = Si + Ni is really a sum of the signal Si plus noise TVj. The noise contributes to 
the auto-correlation such that Mi x Mi = Sf + Nf, but is minimized in a cross -correlation 
Mi x M 2 = (Si + Nx) x (S 2 + N 2 ) = Sf, since Si = S 2 . Another advantage to a cross- 
correlation study in this analysis is the availability of multi-epoch data over a four-year period. In 
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Figure S 3: The cross power spectra of the sum of multi-epoch maps. The cross -correlation 
power spectra of different epoch summed maps with 3.6 /im (top) and 4.5 /mi (bottom) shown 
separately. The average of the summed maps are taken to be the power spectrum. The notation 
(a + b) x (c 4- d) indicates a cross correlation between the average of the a + b and the c + d epochs. 
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such data any time varying signals that are not correlated between epochs cannot contribute to the 
background anisotropy power spectrum. One such possibility is the zodiacal light associated with 
scattered Sunlight off of dust particles. This is due to dynamical dust particles in near-Earth orbits. 
Furthermore with varying Spitzer orbit the lines of sights to the SDWFS Bootes field during the 
four epochs will also be different. Thus, we expect the zodiacal light contamination to have a time 
varying component that is not correlated between epochs. While autocorrelations in single epochs 
may be contaminated by zodiacal light, we expect cross-correlations to reduce the contamination. 
A previous analysis showed that the contamination from zodiacal light spatial fluctuations is at 
least an order of magnitude below the background anisotropy level at 3.6 fim. Thus, zodiacal light 
should not be the dominant systematic effect in the present analysis. The cross-correlations using 
sum maps of epochs 1 to 4 are shown in Fig. $31 

We can check the assumption that Ni x N 2 terms cancel by examining the cross-correlation 
between, say, Mi — M 2 = N± — N 2 and M 3 — M 4 = N 3 — iV 4 . For differences of the same 
region the signal terms cancel and the amplitude of the cross-correlation (Mi — M 2 ) x (M 3 — Ml) 
provides an estimate for the floor level at which noise contributions cancel. In Fig. $4] we show 
these differences, where we find that the cross-correlation power spectra between different epoch 
differences are consistent with zero. The variance of these difference cross-correlations provide a 
part of the error budget associated with noise correlations between different epochs including any 
systematic effects that are not canceled out in difference maps and are correlated between epochs. 
Since these spectra represent the noise floor the absolute value, we add them to the final error 
budget in quadrature with other uncertainties (Fig. SO). 

Even after cross-correlation, the raw spectra are contaminated by several different sources 
and require additional corrections. These issues include resolution damping from the beam, the 
fictitious correlations introduced by the mask and shot noise. All three of these contaminants are 
dealt with as described below. 



4 The Beam Correction 



The first correction that needs to be applied to the raw power spectra is the correction for the beam. 
Realistic detectors have limits to their resolving power which causes a fictitious drop in power 
at high multipoles. For a known beam structure, the resolution limits of the instrument can be 
modeled in harmonic space with a function that encodes the full-width-half-max (FWHM) of the 
telescope. The resulting scale dependent function is known as the beam transfer function, For a 
Gaussian beam it can be computed analytically as 

h = exp(ZVb eam /2) and (3) 

#FWHM ... 

abcam ~ V8E2' 
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where #fwhm encodes the full-width-half-max of the instrument's resolving power. The Spitzer 
team has measured FWH m = 1.9 arcsec. The beam transfer function can also be computed directly 
from the data by measuring the point spread function. Each bright point source in the sky really 
shines a thin point like beam of light that should only illuminate a single pixel of the detector. 
However, due to the finite resolution of the telescope, the source is spread out over many pixels 
and often has a complex shape very different from a Gaussian. Fig. S(6] shows this is true of the 
SDWFS PSF, and for this reason the beam transfer function needs to be calculated directly from 
the PSF. 



In general, the beam transfer function b[ in Fourier or multipolar space is 



/-f-Mpsf 

b f = -IT— (5) 

1 £,-M p oint 



M r 

where C l ps is the power spectrum of M ps f, the observed image of a point source including the 
effects of the telescope, and C l pomt is the power spectrum of a true point source where all the 
light lies in one pixel of the map. Fig. $7] shows the beam calculated using Eq.[5]and the publicly 
available models of the Spitzer-IRAC PSF. The PSFs differ for each epoch and waveband and the 
appropriate beam transfer function for a cross-correlation between maps M\ and M 2 becomes: 



b} x2 = ^Jbjbf (6) 
where bj and bf are the beam transfer functions for maps Mi and M 2 respectively. 



The beam corrected spectra C\ are then computed from the raw power spectrum C\ aw by 
dividing by the beam transfer function bi, 

d = C{ aw /bl (7) 

To measure the uncertainty in the beam transfer function we must understand the uncertainties in 
the PSF. The SDWFS team provides several measured PSFs taken across the Bootes images from 
which the uncertainty in b t can be measured from Eq. [5] as 5b t = 5C™ psl /C z Mpoint with 5C™ psi 
estimated from the variance of the differences from using the various PSF models. 



5 The Mode Coupling Correction. 

Fictitious correlations introduced by the mask must be corrected. When an image is masked, the 
sources are replaced by the value zero in the image. When the power spectrum in computed, these 
zeros in real space make fictitious contributions to the two-dimensional Fourier transform that are 
then added to the final power spectrum. This can be easily seen in Fig. $9] On the top we see an 
unmasked fluctuation pattern for a specific /-mode. After the mask is applied, this fluctuation gets 
broken up into fluctuations of different sizes causing both a diminishing and a reshuffling of power 
in Fourier space. 
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Table S 1 : The final SDWFS power spectrum values 1 2 Ci/2tt for both the 3.6 /iim and 4.5 
fim bands. The quoted error is the 1 a uncertainty of the final power spectrum. 
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A matrix correction method exists^ to model the effects of the mask on the power spectrum 
by using a matrix Mw whose inverse removes the effects of the mask from the measurement by 
matrix multiplication. If Ci is the masked sky power spectrum and C\ be the true power spectrum, 
the relation between the true and masked sky spectrum is 

Q = M IV C V , , (8) 

where Einstein summation notation is being used. Since this relationship is matrix multiplication, 
the masking effects can be removed from the masked power spectrum by simply using a matrix 
multiplication M^, l Ci to recover the true power spectrum C\. In the limit of no /—mode coupling, 
Mui — /sky where / s k y is the fraction of the masked map that is non-zero. 

Calculating the mode coupling matrix Mn> analytically^ is computationally expensive for 
large maps. For this reason, we developed a new way to generate the mode-coupling matrix as 
follows: 



1. For each £ in the power spectrum create many realizations of maps consisting of a pure tone 
where C m = 1 if £ = m and C rn = otherwise (an example case is shown in Fig. $9]). 

2. For each of these trial maps, mask the maps and calculate an observed power spectrum 

C m {£). 

3 . The mode coupling matrix Mi m = (C m {£) ^ is the average of the masked power spectra found 
for the random realizations of model £. The inverse of the mode coupling matrix gives the 
sky power spectrum corrected for the masking effect, Ci = Mj^C m . 



To see how well this this works, consider Fig. $9j The black line shows the exact power 
spectrum from which we drew 100 simulated images. The red points show the observed power 
spectra of the masked sky. The blue points show those same 100 power spectra after correcting with 
the mode coupling matrix described above. This mode-coupling transformation does an excellent 
job recovering the input power spectra of an unmasked sky. As the lower panel of Fig. $9] shows 
using the simplified model based only on the masked sky fraction / S k y is not a good approximation. 



To estimate the effects of cosmic variance on the power spectrum one must make multiple 
Gaussian simulations with the exact power spectrum measured in the data. Without applying a 
mask, these simulations will have a cosmic variance: 

where §Cf y is the cosmic variance for multipole I, SI is the width of that I bin, / s k y is the fraction 
of the total sky covered by the unmasked region, and I in the denominator is the mid-point of the 
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I bin. It should be clear that masks should increase cosmic variance because they further reduce 
the coverage of the unmasked sky. The mode-coupling matrix does not restore errors from cosmic 
variance and therefore, masking and correcting with [Mu>)~ 1 will combine errors from the mode- 
coupling matrix as well as the reduced coverage from the mask and thus have a greater variance 
then cosmic variance alone. 

6 The Map-making Transfer Function 

The last correction that must be accounted for is to correct for the effects of the map-making 
procedure in producing the final mosaic. The mosaicing procedure is often combines many images 
and between dealing with overlap regions as well adding more tiles to a mosaic, the addition or 
subtraction of total power must be quantified and corrected. The best way to uncover the effects 
of the mosaicing procedure on the power spectrum is to take a map with a known power spectrum, 
re-build the map using your mosaicing procedure and compare the output and input power spectra. 
More specifically to build the map-making transfer function we followed the following procedure: 

1. We created a large map similar in size, pixel scale and astrometry to our SDWFS data mosaic, 
but assuming a shape for the power spectrum (we started with pure white noise for which 
Ci = constant). 

2. We break the map into the tiles similar in size, pixel scale and astrometry to our original 
data tiles. We added instrumental noise to the tiles consistent with the noise of the different 
epochs. 

3. We ran these simulated signal plus noise tiles through the same self-calibration mosaicing 
procedure described in Section 1 to produce a new map identical to the original map modulo 
the mosaicing effects. 

4. The map-making transfer function is then 7] = (7° rig /(7 ; mosiaccd ? where C° llg is the known 
power spectrum of the original simulation and Q n0Slaced is the power of the final map. 

5. We repeated Steps 1 to 4 above for different initial power spectra to test if the transfer func- 
tion remains the same or is different. We found that the transfer function is independent of 
the assumption for the input power spectrum shape or the amplitude. 

For our analysis we used this process with the self-calibration algorithm for the mosaic trans- 
fer function. Simulated maps of pure white noise broken into tiles and remosaiced using the 
self-calibration algorithm to determine the transfer functions. An example transfer function for 
SDWFS is given in Fig. S flOl This was obtained by simulating 10 3 independent maps and using 
the mean and the standard deviation of (Tj)i, where i denotes each simulation, to determine the 
best-determined transfer function and its error. 
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Given the transfer function TJ, beam correction hi, the mode-coupling matrix Mu> the final 
power spectrum is estimated as 

Q = M^TlCvjb], (10) 

where Cy is the raw power spectrum after masking out the foreground sources and C\ is the fi- 
nal corrected power spectrum. C\ is the power spectrum that we present in the main paper and 
compared to previous results. 

7 Final Power Spectrum Results 

To compute the power spectrum for the infrared background at 3.6 fim we used all four epochs of 
the SDWFS data as described above. In order to minimize contamination from instrumental noise, 
zodiacal light and other systematics we used the cross-correlation 1/2(E 1 + E 2 ) x l/2(E 3 + E4) 
and the two additional permutations by switching the epochs, where E\..E± correspond to the four 
epochs. Masks were generated to remove the foreground and applied as described above before the 
cross-correlations were taken. We used the same cross-correlation procedure to obtain the 4.5 fim 
power spectrum. After the raw spectra are obtained from the mosaics, they are corrected for the 
beam, mode-coupling, mosaicing (or map-making) effects as described above using Equation [TOl 

Fig. SfTTI shows the results from the Spitzer SDWFS Bootes field (the power spectra values 
are listed in Table [T]). Results from a recent analysis^ are also given for reference. These spectra 
are the final spectra after all corrections have been applied. We note the strong agreement between 
our measurements and the previously published ones. The difference at small angular scales, high 
i, is due to differences in the depth of the mask. It is captured by a difference of the shot-noise 
levels in the point source detection level between SDWFS and deeper SEDS dataPl 

Fig. S fT2l shows the angular cross-power spectrum of near-IR anisotropics measured with 
SDWFS at 3.6 and 4.5 /mi. The left panel shows the cross power spectrum (Cf • 6_4 - 5 ) between the 
two channels, while the right panel shows the correlation coefficient calculated as 

r = Cf^/y/C^Cf- 5 , (11) 

where Cf' 6 and Cf' 5 are the auto power spectra at 3.6 and 4.5 /im, respectively (Fig. SfTTI). The 
correlation coefficient is consistent with unity. The errors are the la overall uncertainty in the 
correlation coefficients found by propagating errors on the cross power spectrum and auto power 
spectra through equation ITTl 

8 Theoretical Interpretation of Near-IR Anisotropics as Spatial Fluctuations of Integrated 
Intrahalo Light 

Our intrahalo light (IHL) model presented in the Letter is described in this Section. In Fig. 1 of 
the main Letter, we also show results from two descriptions related to the near-IR background 
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anisotropics. One involves the faint galaxies that fall below the magnitude cut-off of the masks 
that are applied to measure fluctuations. The contribution from faint galaxies, primarily dwarf has 
been studied in detail with latest information on the faint-end of the galaxy luminosity functions^ 
and we use their results in Fig. 1. A second model involves the z > 6 galaxy contribution. The 
shaded region for z > 6 galaxies in Fig. 1 combines the analytical models 19 with results from 
numerical simulations 20 . The predictions are normalized to the measured luminosity functions of 
galaxies at z > 6 and uses a reionization history that is consistent with the WMAP 7-year optical 
depth to electron scattering^. 



For the interpretation presented in the Letter, we model the IHL intensity angular power 
spectrum using the halo model. The IHL model differs from galaxy clustering models in that we 
assign a profile to the diffuse stars. The standard galaxy clustering models assume a central galaxy 
at the halo center and satellites that are distributed randomly in the halo with a profile that is tracing 
the dark matter distribution. We now add a diffuse extended component in addition to the central 
and satellite galaxies. The halo number density^ as a function of redshift and mass dn(M, z)/dM 
is 

dn p m . dv 



dlnM M J y ' dlnM 



with 



uf(u) = A\ -av 2 [l + {av 2 )~ p } exp[-az/ 2 /2] . (13) 



7T 



Given a galaxy luminosity-halo mass relation, a certain fraction ficL(M) of this luminosity will 
be in the form of IHL. The IHL luminosity-mass relation is then: 



J IHL 



A (M, z) = f IHL (M)L(M, z = 0)(1 + z) a f x (\/(l + z)) , (14) 



where a is the power-law index that accounts for a possible redshift evolution and f\(\/(l + z)) 
is the spectral energy distribution of the IHL (see also discussion below). We model the fraction of 
total luminosity in form of IHL as a power-law in halo mass, 

WMH^U^ . (15) 



The total luminosity as a function of halo-mass L(M, z = 0) at z = is taken to be the 
best-fit relation frorrP^at 2.2 /im 

/ M \ - 72 
L(M, z = 0) = 5.64 -10^( 27 . 1014WM J L Q . (16) 

This is measured for galaxy groups and clusters. At smaller mass scales one no longer has the issue 
of multiple galaxies in a halo and the total luminosity is simply that of the central galaxyS3. We 
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extend it to lower masses using the same power-law slope since analyses of the total luminosity- 
halo mass relation using galaxy-galaxy lensing find the slope continues down to mass scales below 
10 11 MqEsEO. Given that this total luminosity-halo mass relation is measured at 2.2 /im, we scale 
this to other wavelengths with the SED fx, but with the normalization that fx = 1 at 2.2 yum at 
z — 0. We consider the SED of IHL to be consistent with that of old elliptical galaxies comprised 
of old, red stars^. In the main Letter, in Fig. 3, we also consider alternatives for the SED using a 
variety of galaxy SED templates. 

Under these assumptions the angular power spectrum of the IHL flux fluctuations can be 
written as the sum of a 1-halo term, that originates from small-scale fluctuations within individual 
halos 

°l h = 7A2 / dV TTZ I u \ f Mma * dM^^uj RL (k\M)Lj HL ^(M, z) , (17) 
[An) 2 J (1 + zyx (z) JM mm dM 

and a 2-halo term, related to the large scales dark matter fluctuations, and hence to the linear dark 
matter power spectrum Pm(&, z) as 

r \ -i 2 

Cf = -±-(dV 1 [ MmaX dM^^u mL (k\M)b h (M^)L IHL ,x(M,z 
{Any J (1 + z) 2 x [z) [JM min dM 

xP M (k = £/ X (z),z), (18) 

where uihl(&, z\M) is the Fourier transform of the IHL profile in a dark matter halo of mass M at 
redshift z, bh(M) is the linear bias, x( z ) is the comoving radial distance, and dV is the comoving 
volume element dV = x( z ) 2 dx/dz. The redshift integration is performed up to a maximum 
redshift z max = 5. We found that integrating to a higher redshift did not change our results. The 
values of M min and M max in Eqs . [171 and [TBI determine the relative amplitude of the 1-halo and 2- 
halo terms and we let them vary freely. The power spectrum also contains a shot-noise contribution 
from unresolved fluctuations, so that the total power spectrum is 

Ct = C} h + Cj h + Cf . (19) 

We also take C| N to be a free parameter that is varied during our model fitting. 



Since the IHL profile for small mass halos has yet to be determined precisely, we consider 
two model descriptions under the assumption that IHL (i) traces the Navarro-Frenk- White (NFW) 
profile of dark matter halos^l and (ii) falls as r~ 2 with an exponential cut-ofP^ such that p mL oc 
1/r 2 exp(— r/2i? vir ). There is limited information on the light profile from the stacking analysis 
of SDSS galaxies 53 . However, we are unable to use those measurements for the current study as 
we do not have information on how the profile changes with the halo mass. Moreover given the 
limited information both in terms of the angular scale of fluctuations and the large uncertainties we 
find that we are not able to statistically distinguish one IHL profile over another. 
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In Fig. 1 of the main Letter, the best-fit anisotropy power spectrum makes use of the descrip- 
tion involving the NFW profile with 

= 7 — 7 vT^ 7 ^ ' (2 °) 

(cr/r vir )(l + cr/r viT y 

where r vir is the virial radius and c the halo concentration parameter. We define the concentration 
using the result from numerical simulations^ that find 

where M* is the mass scale at which the critical density contrast S c required for spherical collapse is 
equal to the square root of the variance in the initial density field a(M*) = S c . While we make use 
of this particular fitting functional, an alternative fitting function 40 led to the same results. Thus, 
our best-fit parameter values are independent of the assumption on the halo mass-concentration 
relation. 



To analyze the data we conducted a Monte Carlo-Markov Chain (MCMC) analysis using a 
modified version of the publicly available package cosmoMCpSwith convergence diagnostic based 
on the Gelman-Rubin statistic 41 . We fit a total of six free parameters: the minimum and maximum 
masses (M min , M max ), the power-law index of the redshift dependence in the luminosity-mass 
relation a, the amplitude and the power-law index with halo mass of the IHL fraction Af and 
/3, respectively, and the shot noise contribution Cf N . We fix the cosmological parameters to the 
best-fit for the AC DM concordance model from WMAP 7-year analysis^. 



The best-fit parameters for the fit to the data at 3.6 /im from SDWFS are shown in Tabled 
The minimum and maximum halo masses of the halo mass range contributing to measured near- 
IR anisotropics is io 9 03±0 05 M Q and io 1L91±0 05 M , respectively. Our model-fitting suggests a 
small dependence of the IHL fraction on halo mass with a power-law slope of (3 = 0.09 ± 0.01 . We 
found values consistent with the la uncertainties when using the alternate IHL profile described 
above and fitting the model to the 4.5 /im power spectrum measurements. Our results for fmh(M) 
are summarized in Fig. 2 of the main Letter. There we also compare our determination to an 
analytical prediction of the IHL fraction relative to the total luminosity, in the literature^. These 
analytical models exist only down to a halo mass scale of about 5 x 10 10 M and /ihl > 5 x 
10~ 4 . In addition to the power-law behavior it was found also found in previous analytical workPl 
that the IHL fraction is constant at a value of about 0.005 when M < 5 x 10 11 M in some 
scenarios to generate IHL. This flattening behavior could be related to our observation that the 
IHL fraction is not strongly halo mass dependent over the mass ranges that we are probing with 
near-IR background anisotropy power spectrum. More detailed studies are necessary to properly 
understand how our results can be used to understand the merger rate and generation of IHL in low 
mass halos at redshifts of 1 to 4. 



In Fig. SfT3lwe summarize the redshift dependence of the IHL power spectrum by calculating 
dCi/dz as a function of redshift. The contributions peak at a redshift of 3, but has a broad distri- 
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bution ranging from 1 < z < 4. The measured shot-noise level of 120 ± 10 nJy nW m~ 2 sr _1 is 
about a factor of 2 higher than the shot-noise in the deeper SEDS data 12 with a value of ~ 57 nJy 
nW m~ 2 sr™ 1 . 

In Fig. 3 of the main Letter, we summarize the SED of IHL. Here, we make use of a variety 
of stellar templates from B to K-type stars for this purpose. The prediction converges at the longer 
wavelengths due to existing measurements but we find large deviations at 1 /im and shorter wave- 
lengths. A measurement of the background anisotropy at optical wavelengths is clearly desirable 
and could be used to both identify the stars that are primarily contributing to IHL at z ~1 to 4. 
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Figure S 4: The cross power spectra of the difference of multi-epoch maps. The cross- 
correlation power spectra of the difference of multi-epoch maps between epochs 1 to 4 with 3.6 
/im (top) and 4.5 /im (bottom) shown separately. The cross-correlations are consistent with zero 
and the variance between the different cross-correlations provide one part of the final error budget 
associated with the power spectrum measurement. 
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Figure S 5: The noise floor to detect near-IR background anisotropics. The mean and the 
variance (la) of the different-epoch difference maps cross-correlations. We show the absolute 
value of the mean and the error on the mean from the variance of the three independent cross- 
correlations. 
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Figure S 6: The IRAC PSE The 3.6 /im (left) and 4.5 /im (right) IRAC PSFs on a logarithmic 
intensity scale. 



io°r 



10" 



-a 10" 2 r 



10" 3 r 



10" 4 



10' 



10 ; 



10' 



A 3.6/i.m 

x 4.5/zm 



10' 



10 fc 



Figure S 7: The IRAC beam function. The beam transfer functions for the SDWFS Bootes field 
in the multipole space t. We show the functions at 3.6 and 4.5 /im with triangles and crosses, 
respectively. The plotted error bars are the 1 a uncertainties. 
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Figure S 8: The masking effects in a map. The left shows a map of a / = 285 fluctuation. On the 
right is the same map masked by one of the SDWFS masks. Note how what was only a large scale 
fluctuation gets broken up into smaller modes by the mask, contaminating the true power. 



Af 0.0015 ±0.0002 

log(M min /M ) 9.03 ± 0.05 

log(M max /M ) 11.91 ±0.05 

(3 0.094 ±0.005 

a 1.23 ±0.09 

Cf N (nW 2 rrr 4 sr 1 ) (9.8 ± 0.5) x lO" 11 



Table S 2: The best-fit parameter values of the IHL anisotropy power spectrum model to 
the 3.6 /im data using MCMC model fits. The quoted error bars are the \ o uncertainties 
for each of the parameter likelihoods marginalizing over other parameters. 
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Figure S 9: Masking effects on the power spectrum. The top panel shows how well we can 
recover the true spectrum of a masked sky using the mode coupling matrix. The black points show 
the seeded spectra we used to create 100 simulations. The red shows what happens to the spectra 
after masking using the 3.6 fxm mask. The blue points show what happens after correcting with 
the mode-coupling matrix. The bottom panel shows the result if corrected with only the masked 
sky fraction given by / sky . The mask breaks many large modes into smaller modes, so that after 
the f s ky correction the large-scale modes are under-represented and the small-scale modes are over 
represented. This illustrates the need to use the^full coupling matrix to correct for the mask. The 
plotted error bars are the 1 a uncertainties. 
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Figure S 10: The map-making transfer function. Spitzer-IRAC mosaic transfer function based 
on the self-calibration algorithm used to make data maps for this study. We show the transfer 
function at 3.6 /im. The 4.5 fxm transfer function is similar to the result at 3.6 /zm. The plotted 
error bars are the 1 a uncertainties. 



22 



10°f 



10" 1 r 



10" 



P 10 _4 r 



10" 



C, = 6.2x10 
C, SN = 4.4x10 




/ .* 
■ *•< 



• SDWFS 3.6/i.m 

x Kashlinsky et al. 201 2 

■ ■ ■ ■ I 1 1 I I u 



10' 



10' 



10' 



10' 



in 



CN 



10°f 



10" 1 r 



0" 



1 0" 3 t 



P 10- 4 r 



10" 



C, SN = 3.9x10-" 
C, SN = 2.0x10"" 




• SDWFS 4.5/i.m 

x Kashlinsky et al. 2012 

■ ■ ■ ■ i i i i i i_ 



10' 



10' 



10' 



10' 



Figure S 1 1 : The angular power spectrum of near-IR anisotropies. The angular power spectrum 
of near-IR anisotropies measured with SDWFS at 3.6 and 4.5 /mi. The 1 a error bars include all 
uncertainties we have discussed in the Supplement and the measurements are beam corrected. 
We also compare our measurements to existing results^ where we find a general agreement on 
clustering. The large-£ difference between the two datasets reflect the depth of the point source 
identification and removal in the mask. 
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Figure S 12: The angular cross-power spectrum of near-IR anisotropies. The angular cross- 
power spectrum of near-IR anisotropies measured with SDWFS at 3.6 and 4.5 /im. The 1 a error 
bars include all uncertainties we have discussed in the Supplement and the measurements are beam 
corrected. The upper panel shows the cross power spectrum (Cf 6 ~ 4 ' 5 ), while the lower panel shows 

the correlation coefficient calculated as r = Cf 6 ~ 4 ' 5 /wCf- 6 C; 4 - 5 , where Cf 16 and Cf' 5 are the auto 
power spectra at 3.6 and 4.5 /im, respectively (Fig. S fTlj) 
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Figure S 13: The redshift dependence of the IHL anisotropy power spectrum. dCi/dz as a 

function of redshift for t = 3 x 10 3 and 10 4 . The majority of near-IR anisotropics originate from 

1 < z < 4. 
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